So far you will be familiar with the binomial model and the Black Scholes model for option pricing
Fundamentally as a result of the Cameron Martin Girsanov Theorem and the Central Limit Theorem both these methods are effectively the same thing
As the number of steps tends to $\infty$ the binomial model tends to the Black Scholes model
The Monte Carlo methods are also fundamentally the same thing - but with the Monte Carlo methods we can start from a binomial perspective or from a Black Scholes perspective and gain insight into how all the methods are effectively equivalent
Additionally with the Monte Carlo method it is much easier to extend the method to deal with more advanced modelling of equity prices with features such as skew and stochastic volatility
Remember when you first looked at the binomial model
You started with an uptick and a downtick
Then do some algebra and the formula that emerges is $f_{n} = e^{-r\Delta t} \times E_\mathbb{Q} \left(f_{n+1}\right)$, that is effectively an expected value over some 'other' probability measure '$\mathbb{Q}$'
We start with a very simple model of the price evolution of our underlying stock
The stock has price $S_0$ at time 0 and the model advances in units of time $\Delta t$
The term of the option is $T$
We assume the volatility of the stock is $\sigma$, which leads to a natural choice of $u$ and $d$ of:
$u=e^{\sigma \sqrt{\Delta t}}$ and $d=\frac{1}{u}$
and $q=\frac{e^{r \Delta t}-d}{u-d}$ - from arbitrage free construction
The value at each node is then worked out from the subsequent nodes with the formula: $f_{n} = e^{-r\Delta t} \times E_\mathbb{Q} \left(f_{n+1}\right)$
These formulas are derived by considering a hedging portfolio of stocks and cash
The important thing about $u$ and $d$ is the gap between them as the Girsanov Theorem can sort out the drift but not the volatility
The $\sqrt{\Delta t}$ is to make the unit time volatility $\sigma$
So simply by changing the probabilities of the uptick and the downtick we can express the market price of an option (or any other financial asset) as an expected value
This extends through the Cameron Martin Girsanov Theorem to the full Black Scholes Model.
An illustration of this can be seen in this graphical simulation
Much of the mathematics of the binomial method has to be reproduced to use the Monte Carlo method
although Monte Carlo methods are often simpler than other methods we CANNOT
All we can do in the Monte Carlo method is simplify the algorithm by being able to pick the up tick and the down tick randomly rather than having to write the code to go through every single permutation
For given values of $S_0$, $K$, $T$, $\Delta T$ and $\sigma$ can you write a Monte Carlo version of the binomial model
Can you explain how this is simpler to code than the normal binomial model
There is just one path for each iteration
What other advantages does it have?
You do not run the risk of running out of memory especially for those options which are not recombining
Fundamentally the CMG is about $\frac{d \mathbb{Q}}{d \mathbb{P}}$ but the visible impact of this is to appear to shift the distribution from the expected return on equity to the risk free rate as we have seen in the above demonstration
Can you change your model to reflect this directly in the VBA
Hint: Preserve the ratio of $u:d$ and solve $q=0.5$
Solution:
$u=C \times e^{\sigma \sqrt{\Delta t}}$ and $d=C \times e^{-\sigma \sqrt{\Delta t}}$
$q=\frac{e^{r \Delta t}-d}{u-d} = 0.5$
$q=\frac{e^{r \Delta t}-C \times e^{-\sigma \sqrt{\Delta t}}}{C \times e^{\sigma \sqrt{\Delta t}}-C \times e^{-\sigma \sqrt{\Delta t}}} = 0.5$
$C=\frac{e^{r \Delta t}}{e^{-\sigma \sqrt{\Delta t}} + 0.5 \times e^{\sigma \sqrt{\Delta t}} - 0.5 \times e^{-\sigma \sqrt{\Delta t}}}$
$C=\frac{2 \times e^{r \Delta t}}{e^{\sigma \sqrt{\Delta t}} +e^{-\sigma \sqrt{\Delta t}}}$
It is a necessary limitation of the binomial model to only have discrete state spaces.
This clearly means an uptick and a downtick in the case of the binomial model but a basic extension is the trinomial model which has three ticks from each state
However for a Monte Carlo approach we are not limited to discrete state spaces.
Why
because we do not store each state
Therefore it is possible (and sensible) to extend the model by allowing the return at each point to be generated fully from a lognormal distribution
Extend your model to allow for full lognormal steps
How many steps is it now sensible to have to get to an accurate result (You may assume sufficient Monte Carlo runs).
ONE as we only needed multiple steps under the binomial model to allow time for the central limit theorem to get us to a lognormal distribution
For simple options with lognormal steps ONE step is the Monte Carlo equivalent of the Black Scholes model and will converge to the Black Scholes result as the number of runs becomes large
Test your model using just one step
Can you reconcile the mathematics of what we have done to the proof of the Black Scholes formula
The Monte Carlo now is simply the final stage of the Black Scholes proof - that is the integration at the end once we have established that the arbitrage free value is genuinely the expected value over $\mathbb{Q}$
If you have followed the steps above you will probably expect to have replicated Black Scholes results with a one-step lognormal model.
However (unless you have thought ahead to this section) it is highly likely you will have made a mistake and your answers will be too high
The mistake we have made is to assume that the solution of the integral equation: $\frac{dS_t}{S_t}=r dt + \sigma dW_t$ is $S_t=S_0 e^{rt+\sigma W_t}$
When in fact it is $S_t=S_0 e^{(r-0.5 \sigma^2)t+\sigma W_t}$
$\frac{dS_t}{S_t}=r dt + \sigma dW_t \implies dS_t=r S_t dt + \sigma S_t dW_t$
Now look at Ito's Lemma:
$df(t, X_t) = \left[\frac{df}{dt} + \mu (t, X_t) \frac{df}{dx} + \frac{1}{2} \sigma^2 (t, X_t) \frac{d^2f}{dx^2} \right] dt + \sigma(t, X_t) \frac{df}{dx} dW_t$
Where $dX_t=\mu(t,X_t)dt + \sigma(t, X_t) dW_t$
We set this up as follows:
Rewrite with $S_t$ instead of $X_t$
Then our $\mu$ is $r S_t$ and our $\sigma$ is $\sigma S_t$
The we use the function $f(t, S_t)=ln(S_t)$ (there is no way of knowing this other than it works)
so $df(t, S_t) = \left[\frac{df}{dt} + \mu (t, S_t) \frac{df}{ds} + \frac{1}{2} \sigma^2 (t, S_t) \frac{d^2f}{ds^2} \right] dt + \sigma(t, S_t) \frac{df}{ds} dW_t$
then $df(t, S_t) = \left[\frac{df}{dt} +r S_t \frac{df}{ds} + \frac{1}{2} \sigma^2 S_t^2 \frac{d^2f}{ds^2} \right] dt + \sigma S_t \frac{df}{ds} dW_t$
Now we need to calculate our derivatives:
$f(t, S_t)=ln(S_t) \implies \frac{df}{dt}=0, \frac{df}{ds}=\frac{1}{S_t}, \frac{d^2f}{ds^2}=\frac{-1}{S_t^2}$
So putting it all together:
$d ln(S_t)= \left[0 +r S_t \frac{1}{S_t} + \frac{1}{2} \sigma^2 S_t^2 \frac{-1}{S_t^2} \right] dt + \sigma S_t \frac{1}{S_t} dW_t$
Simplify
$d ln(S_t)= \left[r - \frac{1}{2} \sigma^2 \right] dt + \sigma dW_t$
Integrate:
$\displaystyle\int d ln(S_t)= \displaystyle\int (r - \frac{1}{2} \sigma^2) dt +\displaystyle\int \sigma dW_t$
$ln(S_t)= (r - \frac{1}{2} \sigma^2)t + \sigma W_t + ln S_0$
$\therefore S_t=S_0 e^{(r-\frac{\sigma^2}{2})t + \sigma W_t}$
We can now see below the code for the three different methods
Option Explicit Function max(x As Double, y As Double) As Double If x > y Then max = x Else max = y End If End Function Function CRR(S0, K, r, T, sigma As Double, n As Integer) As Double Dim S() As Double ReDim S(0 To n, 0 To 2 ^ n - 1) Dim f() As Double ReDim f(0 To n, 0 To 2 ^ n - 1) Dim delta_t, u, d As Double Dim j As Double Dim q As Double Dim node_x, node_y As Integer delta_t = T / n u = Exp(sigma * delta_t ^ 0.5) d = Exp(-sigma * delta_t ^ 0.5) q = (Exp(r * delta_t) - d) / (u - d) S(0, 0) = S0 For node_x = 1 To n For node_y = 0 To 2 ^ node_x - 1 If node_y Mod 2 = 0 Then S(node_x, node_y) = S(node_x - 1, node_y / 2) * u Else S(node_x, node_y) = S(node_x - 1, (node_y - 1) / 2) * d End If Next node_y Next node_x For node_y = 0 To 2 ^ n - 1 f(n, node_y) = max(S(n, node_y) - K, 0) Next node_y For node_x = n - 1 To 0 Step -1 For node_y = 0 To 2 ^ node_x - 1 f(node_x, node_y) = Exp(-r * delta_t) * (q * f(node_x + 1, node_y * 2) + (1 - q) * f(node_x + 1, node_y * 2 + 1)) Next node_y Next node_x CRR = f(0, 0) End Function Function Monte_Carlo_Binomial(S0, K, r, T, sigma As Double, n, runs As Integer) As Double Dim delta_t, u, d As Double Dim j As Long Dim q As Double Dim S As Double Dim sum As Double Dim payoff As Double Dim run As Long delta_t = T / n u = Exp(sigma * delta_t ^ 0.5) d = Exp(-sigma * delta_t ^ 0.5) q = (Exp(r * delta_t) - d) / (u - d) sum = 0 For run = 1 To runs S = S0 For j = 1 To n If (Rnd() < q) Then S = S * u Else S = S * d End If Next j payoff = max(S - K, 0) sum = sum + payoff Next run Monte_Carlo_Binomial = Exp(-r * T) * sum / runs End Function Function Monte_Carlo_LogNormal(S0, K, r, T, sigma As Double, n, runs As Integer) As Double Dim delta_t, u, d As Double Dim j As Long Dim q As Double Dim S As Double Dim jump As Double Dim sum As Double Dim payoff As Double Dim run As Long delta_t = T / n sum = 0 For run = 1 To runs S = S0 For j = 1 To n jump = Exp(Application.NormSInv(Rnd()) * sigma * Sqr(delta_t) + (r-0.5*sigma^2) * delta_t) S = S * jump Next j payoff = max(S - K, 0) sum = sum + payoff Next run Monte_Carlo_LogNormal = Exp(-r * T) * sum / runs End Function